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Abstract. A frequently faced task in experimental physics is to measure the probability 
distribution of some quantity. Often this quantity to be measured is smeared by a non- 
ideal detector response or by some physical process. The procedure of removing this smearing 
effect from the measured distribution is called unfolding, and is a delicate problem in signal 
processing, due to the well-known numerical ill behavior of this task. Various methods were 
invented which, given some assumptions on the initial probability distribution, try to regularize 
the unfolding problem. Most of these methods definitely introduce bias into the estimate of 
the initial probability distribution. We propose a linear iterative method (motivated by the 
Neumann series / Landweber iteration known in functional analysis), which has the advantage 
that no assumptions on the initial probability distribution is needed, and the only regularization 
parameter is the stopping order of the iteration, which can be used to choose the best compromise 
between the introduced bias and the propagated statistical and systematic errors. The method 
is consistent: "binwise" convergence to the initial probability distribution is proved in absence 
of measurement errors under a quite general condition on the response function. This condition 
holds for practical applications such as convolutions, calorimeter response functions, momentum 
reconstruction response functions based on tracking in magnetic field etc. In presence of 
measurement errors, explicit formulae for the propagation of the three important error terms 
is provided: bias error (distance from the unknown to-be-reconstructed initial distribution at a 
finite iteration order), statistical error, and systematic error. A trade-off between these three 
error terms can be used to define an optimal iteration stopping criterion, and the errors can 
be estimated there. We provide a numerical C library for the implementation of the method, 
which incorporates automatic statistical error propagation as well. The proposed method is also 
discussed in the context of other known approaches. 



1. Introduction 

In data analysis one commonly faces the problem that the probability density function (pdf) 
of a given physical quantity of interest is to be measured, but some random physical process, 
such as the intrinsic behavior of the measurement apparatus, smears it. The reconstruction 
of the pertinent pdf based on the measured smeared pdf and on the response function of the 
measurement procedure is called unfolding. To be specific, let us have the original unknown 
pdf X (-). /(x) of the undistorted physical quantity which we need to reconstruct, and assume 
that the actual measured pdf can be expressed of the form y i— g{y) = f p{y\x) f{x) dx, where 
{y,x) I—)- p{y\x) describes the smearing effect in a probabilistic manner!^ Then, it is said that 

^ All pdfs are understood to be real valued non-negative Lebesgue integrable functions over some finite 
dimensional real vector space X. 



the pdf g is the pdf / folded with the response function plf| Our mathematical task is to solve 
the above linear integral equation in order to obtain /, given g and p. This problem is known 
not to be a simple numerical task (ill-posed problem), and several methods are used by the data 
analysis communities in order to regularize the problem (for an overview on the most popular 
approaches, we refer to [HE])- 

Let us denote by Ap the pertinent folding operator, which acts like {Apf) (y) = 
/ p{y\x) f{x) dx on a function / at a point yH Given the measured pdf g = Apf , the problem 
of unfolding can then be formalized as follows: the pdf / = A~'^{g) is to be determined or 
approximated. The mathematical cause of the numerical ill-posedness of this unfolding problem 
can then be put forward as: the inverse A~^ of a generic folding operator can be shown not to be 
continuous despite the forward folding operator Ap always being continuou^ (this phenomenon 
is discussed in detail e.g. in [9]). The non-continuity of the inverse folding operator A~^ may 
be also reformulated in a less abstract manner: initially distant functions can be mapped close 
by the folding operator Ap, as illustrated in Figure [TJ I.e. one can lose discriminating power 
between pdfs upon a folding. 




Ap(fj) close ApCf^) 



Figure 1. Illustration of the non-continuity of the inverse of a folding operator Ap: two 
distant functions /i and /2 may be mapped close by the folding - distance of functions are 
here understood as probabilistic distance, i.e. in the L^{X) function norm. 

A further aspect of the numerical ill-posedness of the unfolding problem is that in practice 
the folded pdf g is often obtained via statistical measurements (e.g. histograming), and therefore 
is contaminated by statistical errors. I.e. in reality g = Apf + e holds instead of the idealized 
equation g = Apf, where e(x) is a random variable for each point x (or for each histogram bin - 
in the language of histograms) . Thus, when estimating the unfolded pdf as Ap^{g) = f + Ap^{e), 
the contribution of the second term is not guaranteed to remain small due to the non-continuity 
of the inverse folding operator Ap^ even when e is initially known to be small. On top of this, the 
statistical error term e may contain modes not within the image of the folding operator Ap, on 
which the evaluation of the inverse operator Ap^ is not meaningful if the problem is not initially 
discrete. These effects are demonstrated in Figure [21 which shows that simple inversion of the 
discretized folding operator on the measured pdf gives unphysical numerical result: a result very 
different from the initial pdf, having large negative and positive alternating amplitudes. 

In order to regularize the numerical ill-posedness of the unfolding problem, various methods 
are used. These methods can be divided into three large classes. 

Whenever the response function p is translation invariant in the sense that for all x,y,z € X one has 
p{y + z\x) — p{y\x — z), the folding is specially called convolution, and in that case p may be expressed by 
a single pdf: p{y\x) — p{y — x\0). 

^ To be precise, Ap is a L^{X) — >■ L^{X) continuous linear operator, where L^{X) denotes the normed space 
of complex valued integrable functions over the vector space X. The response function p is assumed to be 
p{-\x) e L^{X) for all x £ X. 
* Continuous in the L^{X) — >■ L^{X) sense. 
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Figure 2. (Color online) Demonstration of the numerical ill-posedness of the unfolding 
problem: a Cauchy distribution is convolved with a Gauss distribution with Monte Carlo 
method to generate the measured distribution contaminated with statistical errors. Clearly, 
the unfolded pdf, obtained by simple numerical inversion of the discretized folding operator on 
the measured pdf gives physically unreasonable numerical result: large alternating positive / 
negative amplitude pdf values. 



(i) Using a parametric Ansatz for /, and fit parameters, so that Apf gets close to g. This 
method can be slightly insensitive to the details of the true / (as illustrated in Figured]), 
and of course can introduce strong systematic bias on the result if the parametric Ansatz 
does not hold in an exact manner of the form that was assumed. Such methods are used in 
general for inclusive particle identification by specific ionization (see e.g. |10j). 

(ii) Bin- by-bin fitting of the bin values of the histogramed /, so that Apf gets close to g. This is 
basically equivalent to the naive inversion of the discretized folding operator, and therefore 
produces similar oscillatory results, except when an artificial penalty function is added to 
the in order to suppress large local gradients. In that case, the method can provide 
meaningful answers, but the introduced systematic bias is difficult to quantify. Similarly to 
the parametric Ansatz method, the fit can be slightly insensitive to the details of the true 
/. Most popular methods, such as SVD method |3j, are based on this idea. 

(iii) The iterative method of convergent weights (also known as iterative Bayesian unfolding) 
of Kondor-Miilthei-Schorr-d'Agostini [H [5l El El [8]. This method is, as opposed to the 
previously mentioned methods, is non-linear. On the other hand, by construction it 
preserves positivity and integral of the initial pdf, and therefore maps a pdf exactly into 
a pdf, which does not hold for linear methods, thus, this approach is quite favorable for 
statistical applications. Regularization is achieved solely by stopping the iteration at a finite 
order. However, there is no known proof yet if the iterated pdfs converg^ to the initial 
to-be-reconstructed pdf in a non-discrete scenario, even in the absence of measurement 
errors [6]. Also propagation of statistical and systematic errors of the measured pdf to the 
unfolded pdf has not been investigated, and consequently no generally applicable iteration 
stopping condition is known. 

In a previous paper [9] we proposed a linear iterative unfolding method, for which under certain 
conditions convergence to the initial pdf was proved analytically for some unfolding problems 



^ In case of an iterative unfolding method it is an absolute must to show that the sequence of iterated unfolded 
pdfs converge to the initial one, in the absence of measurement errors (consistency of the method). 



in probability theory (such as convolutions), and due to the linearity of the method, exact 
propagation of statistical errors of the measured (folded) pdf to the unfolded pdf was possible. 
In this paper we propose an improved version of that algorithm, which could be proved to 
be convergent in quite general cases for unfolding problems in a probability theory settinglf] 
The key equality of the convergence proof leads to explicit error propagation formulae for the 
three important error terms: for the bias error (distance from the true unfolded pdf), for the 
propagated statistical error, and most notably for the propagated systematic error, which is 
of great importance in reporting experimental results. An implementation of the algorithm is 
written as C library, along with application examples [13]. The implementation also incorporates 
automatic statistical error propagation. 

The paper is organized as follows: in Section [2] the algorithm and its convergence theorem 
shall be formulated. Section O is devoted to the corresponding error propagation formulae which 
help to formulate an optimal stopping criterion and error estimates therein, while in Section [3] 
we demonstrate our method on examples. 



2. A linear iterative unfolding algorithm 

We provide now a linear iterative solution for a probability theory unfolding problem of the 
form g = Apf, where / is the initial (unknown) pdf, g is the folded (measured) pdf, and p 
is the response function. Given the response function p, one can also define along with the 
folding operator Ap the transpose folding operator by swapping the variables of the response 
function!^ Then, one can attempt to approximate the true unfolded pdf / in the following way: 
define the function sequence by setting the normalization factor 

Kp = max / / p{y\z) p{y\x) dy dz (1) 



and then taking the 

/o = K^'A^pg, 
fN+i = fN+{fo-K;'A'j;ApfN) (2) 

iteration formula. We provide a convergence result on this iterative approximation below in 
absence of measurement errors on g (which is necessary for the consistency of the method). 

Theorem 1. (Convergence) The function sequence N resulting from the above iteration 

scheme converges to the closest possible function to the true unfolded pdf f in the average over 
any compact region, whenever the normalization factor Kp is finite. I.e. for all compact sets 
S <Z X one has 

l^'^oo VoluL(5) Is " " = 

Here, PKer(Ap) denotes the orthogonal projection operator to the kernel set of Ap, and thus 
Pkcv{Ap) = holds automatically whenever Ap is invertible. In addition, the convergence shall 
also hold in the space of square-integrable functions, i. e. one has also 



lim 



f - PKer(Ap)f - fN (x) dx = 0. (4) 



^ The detailed mathematical proof of convergence shall be published elsewhere: [TJ. The proposed iteration 
scheme was motivated by the so called Neumann series and Landweber iteration [12] known in functional analysis, 
but the convergence of neither iterations hold, unfortunately, in a probability theory setting in their original form, 
as one can prove. Our improved iterative algorithm, however, is specially developed to be convergent for unfolding 
problems in probability theory. 

^ The transpose folding operator is defined by {A'^ f){y) — J p{x\y) f{x) dx for all functions / and points y. Note, 
that this simply translates to matrix transposition whenever the folding is discretized. 



Proof The proof is based on Riesz-Thorin theorem and on the spectral representation of 
positive operators in the space of complex square-integrable functions over X (to be published 
in a more mathematically specialized journal: [11]). 

The following observations help to shed some further light on properties of the proposed 
unfolding algorithm. 

(i) In case the pdfs are modeled with histograming, the setwise convergence of pdfs means 
binwise convergence of histograms, i.e. the probability of each histogram bin is restored in 
the limit of infinite iterations. 

(ii) When the inverse of Ap exists, the original pdf / is completely restored. Whenever the 
pertinent inverse does not exist, still the maximum possible information about / is restored, 
namely the function / — i-ker(Ap)/- 

(iii) Whenever Ap is a convolution, then Kp = 1 holds automatically, i.e. Kp < oo is satisfied. 

(iv) The convergence condition Kp < oo holds provably for a wide class of practically relevant 
response functions, such as energy response function of calorimeters, momentum response 
function of track reconstruction in magnetic field etc. 

(v) The iteration scheme of the theorem is motivated by the Neumann series known in functional 
analysis. A similar iterative solution, also referred to as Landweber iteration [T2], is known 
in the theory of Fredholm operators. In probability theory unfolding problems, however, 
the necessary convergence criteria for Neumann series or for Landweber iteration do not 
hold in their original form. 

(vi) The proposed iterative unfolding algorithm does not necessarily need an initial binning 
of pdfs. It may be implemented as well by different density estimators than histograms. 
However, when the pdfs are modeled by histograms, one may recognize that the binning and 
truncation of histograming domain can also be considered as folding operator. Therefore, 
the histogram binning and truncation effect may be included in the response function p, 
and then the effect of histograming can be unfolded (to the maximum possible extent) as 
well. If one wants to numerically implement this, the initial pdf / must be assumed to 
better approach the continuum pdf, i.e. must be assumed to be an unknown histogram 
over a larger domain with finer granulation than the folded one. The schematic of such 
possible rebinning trick is illustrated in Figure [3l 




Figure 3. (Color online) Illustration the rebinning trick for unfolding the histogram binning 
and domain truncation as well (to the maximum possible extent) along with the smearing effect 
of the response function p. For this, the implementation of the folding operator Ap must map 
histograms over a larger domain and with finer graining to histograms with the binning scheme 
of the measured (folded) pdf. 



3. Bias, statistical and systematic errors of the unfolded distribution 

In real measurements, the folded pdf g also admits statistical and systematic errors, and the 
propagation of these terms into the unfolded pdf is necessary to quantify at each finite iteration 
step. The key equality of the proof of Theorem 1 leads to explicit error propagation formulae 
for bias error (distance from the true unfolded pdf), statistical error, and systematic error. First 
we present our result about bias error. 

Theorem 2. (Bias error) Take the iterative solution for the unfolding problem as in Section\^ 
Then, if the normalization factor Kp is finite, the distance of an N-th iterate f^ from the closest 
possible function to the true unfolded pdf f in the average over a compact region has the following 
upper bound: for any compact set S C X one has 



1 



Volume(S') Js 



/ (/ - ^Ker(^p)/ - In) {x) dx 

<j s 



< 



1 



VVolume(5) + 2 



|/|2(x)dx, (5) 



i. e. the residual deviation of /tv from the limit averaged to the set S decreases as . ^ 

Y Volumc(S) 

N+2- 



and 



The above result, translated to the language of histograms means that the bin-by-bin average 
deviation from the true unfolded pdf is bound by the right hand side of the inequality in 
Theorem [JTl where Volume(S') is the histogram bin volume, A'^ is the iteration order, and 



/ |/p(x)dx is an unknown coefficient depending on the true unfolded pdf /. This unknown 
coefficient, however, may be substituted by the calculable expression J f \ f]^\'^{x) dx for large 



N . It is seen that the bias error tends to zero with increasing iteration order N . 

In practical applications, the pdfs are often measured by statistical methods (e.g. 
histograming). In that case, the value of the folded pdf g in each histogram bin admits a 
statistical error. The below theorem states an exact formula for the propagation of this error 
into the unfolded pdf. 

Theorem 3. (Statistical error) Take the iterative solution for the unfolding problem as in 
Section Let C be the covariance matrix of the measured pdf g, where g is assumed to be 
of the form of a histogram. Since a covariance matrix C is positive definite, it is always possible 
to decompose it - not uniquely - in the form C = E for some matrix E, (•)-^ being the 
matrix transpose. (Whenever C is diagonal, construction of such an E is just trivial.) Then, 
the following iteration calculates the statistical error propagation: 



Eq 

En+1 



k~^aIe 



En + (^0 - K^^A^ApEN 



(6) 



where in each step the covariance matrix of fjsf shall be Cn = Ej^ EJ^ . 



Due to the linearity of the method, the contribution of the propagated statistical error term 
is exactly calculable by means of the above formulae, if it is known for the measured pdf g. This 
error term increases with increasing iteration order N. The statistical error of a given histogram 
bin of the A^-th iterate fjy is nothing but the square-root of the corresponding diagonal element 
of C7^. 

Whenever the folded pdf g is a result of an experiment, it may admit a systematic error 6g. 
Also the systematic error 6p of the response function p may give a non-zero contribution to it: 
Aspf. The effect of this initial systematic error on the unfolded pdf is quantified by the following 
theorem. 



Theorem 4. (Systematic error) Take the iterative solution for the unfolding problem as in 
Section 0. Assume that 6g is the systematic error of g (possibly including contribution from 
systematic error of the response function). Then the systematic error for the N-th iterate /at 
averaged over a compact region has the following upper bound: for any compact set S C X 



1 



Volume(S') 



SfNix) dx 



< 



1 



^Volume(S) 



[^{N + 2) + 7) \Kp^Aj6g\^ (x) dx, 



(7) 



i.e. the systematic error of fj^ averaged to the set S decreases as 



and increases as 



^/Volume(5) 

^{N + 2) + 7 ~ 1 + ln(A^ + !)• (^ being the digamma function, 7 being Euler's constant.) 

The above resuh, translated to the language of histograms means that the bin-by-bin average 
systematic error of the A^-th iterate /at is bound by the formula in the right hand side of the 
inequality in TheoremHl where Volume(S) is the histogram bin volume, N is the iteration order, 



and the coefficient J J 



Kp^ATjg 



{x) dx is calculable knowing the bin-by-bin systematic errors 
y, this contribution increases logarithmically with increasing 



5g of the measured pdf g. Clear! 
iteration order A^. 

As the bias error decreases, while the statistical and systematic error of the A^-th iterate /at 
increases with the iteration order N , a trade-off between these error terms provides an optimal 
cutoff criteriorH in the iteration order A^, and error estimates therein. Consequently the true 
unfolded pdf / can be approximated optimally and the error of this approximation can be put 
under full control. Thus, the regularization of the numerically ill-posed unfolding problem is 
achieved, in case of the proposed approach, solely by using an iterative approximation and 
choosing an optimal iteration stop order taking into account the convergent terms (bias error) 
and the divergent terms (statistical and systematic errors). 



4. Examples 

In this section we give two examples to demonstrate our method. In the first example, we take 
a Cauchy distribution, and convolve it with a Gaussian distribution with Monte Carlo method. 
The folded pdf is determined by histograming the sum of the Cauchy and Gaussian distribution 
random numbers, i.e. the measured pdf shall admit Poissonian statistical errors. In this example, 
a relatively modest statistics of 5000 entries was taken to be able to judge the method in the 
low statistics limit. The results is shown in Figure HI It is seen that the original Cauchy pdf is 
restored, modulo the fluctuations arising from the propagated statistical errors ~ these are seen 
as "shoulders" of the unfolded pdf, the amplitude of which decrease with increased statistics. 
The iteration was stopped when the integral of the statistical error term reached about 5% level. 

In the second Monte Carlo example, we generate the energy distribution of transversely 
emitted hadrons in 7GeV p+p collisions [14J, and we assume that this particle spectrum was 
measured by the CMS-HCAL calorimeter [15j. The unfolded spectrum, along with the true and 
measured distribution is shown in Figure [5j 



5. Concluding remarks 

We proposed a linear iterative spectrum unfolding method for application in data analysis. 
Convergence to the true unfolded pdf is proved under a quite general condition [11] in absence 
of measurement errors, and error propagation formulae are derived for bias error, statistical 
error, and systematic error in the presence of measurement errors. The method is demonstrated 
on physical examples. A numerical library in C is provided with the implementation of the 
method [T3]. The algorithm could be included in the ROOUnfold package [16] in the future. 

* E.g. one can take the sum of the three error terms, and stop the iteration when it reaches a minimum. 
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Figure 4. (Color online) Test example with unfolding a Cauchy distribution convolved with 
a Gaussian distribution. Iteration was stopped when the integral of the statistical error term 
reached about 5%. 
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Figure 5. (Color online) A physical example with unfolding energy distribution of charged 
hadrons measured with hadronic calorimeter. Iteration was stopped when the integral of the 
statistical error term reached about 2.3%. 

Acknow^ledgments 

I would like to thank prof. Giinter Zech for valuable discussions. 
References 

[1] Cowan G 2002 Proceedings of Conference on Advanced Statistical Techniques m Particle Physics (18-22 

March 2002, Durham, United Kingdom) (IPPP/02/39, Durham) p 248 
[2] Blobel V 2008 Unfolding for HEP Experiments {Talk at DESY Computing Seminar) 

http : //www . desy . de/~blobel/DESYcompsem08 . pdf 
[3] Hoecker A and Kartvelishvili V 1996 Nucl. Instr. Meth. A 372 469 
[4] Shepp L A and Vardi Y 1982 IEEE Trans. Med. Imag. 1 113 
[5] Kondor A 1983 Nucl. Instr. Meth. 216 177 
[6] Miilthei H N and Schorr B 1987 Mat. Meth. Appl. Sci. 9 137 
[7] Miihhei H N and Schorr B 1987 Nucl. Instr. Meth. A 257 371 



[8] D'Agostini G 1995 Nucl. Instr. Meth. A 362 487 

[9] Laszlo A 2006 J. Phys. A 39 13621 
[10] Laszlo A 2008 Phys. Rev. C 77 034906 

[11] Laszlo A Convergence results on a linear iterative spectrum unfolding method {in preparation) 

[12] Landweber L 1951 Am. J. Math. 73 615 

[13] Laszlo A 2011 The libunfold package {Source code) 

http : //www. rmki .kf ki .hu/~laszloa/dowiiloads/libunf old. tair .gz 
[14] Khachatryan V et al 2010 Phys. Rev. Lett. 105 022002 
[15] Yazgan E et al 2009 J. Phys. Gonf. Ser. 160 012056 
[16] Adye T et al The ROOUnfold package 

http : //hepunx.rl . ac .iik/~adye/software/vmfold/RooUnf old. html 



